%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Program: makePlots.m
% Purpose: make graphics for ALO rust belt project
% By:      David Lagakos
% Date:    05-28-21
% Note:    Only plots vs data when nPeriods = 51
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nPeriods == 51
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% Plot model vs data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%% Figure 1: 2x2 model and data
    figure;
    hold on;
    
    %%% Rust Belt employment share
    subplot(2,2,1)
    RB_labor_share = shiftdim(sum(l_t(1,1:nRBgoods,:)+l_t(2,1:nRBgoods,:)))';
    plot(year',RB_labor_share,'r-',year',mfg_emp_share_RB_data','k:','LineWidth',2)
    axis([1950 2000 0.30 0.55])
    title('Rust Belt Employment Share')
    legend('Model','Data','Location','southwest')
    
    %%% Rust Belt wage premium
    subplot(2,2,2)
    wage_prem_a = shiftdim(sum(rent_t(1,1:nRBgoods,:)));
    wage_prem_b = shiftdim(-kappa*sum(strike_t(1,:,:)./nRBgoods));
    wage_premium_RB = (wage_prem_a + wage_prem_b)./RB_labor_share';
    plot(year',wage_premium_RB,'m--',year',rb_mfg_simple_ratio_data,'k:',year',rb_mfg_dummy_data,'k:','LineWidth',2)
    title('Rust Belt Wage Premium')
    axis([1950 2000 0.00 0.20])
    legend('Model','Data','Location','southwest')
    
    %%% Import shares, model and data
    subplot(2,2,3)
    py_t = p_t.*y_t;
    vtau = tau0.*delta_tau.^(0:1:nPeriods-1);
    RB_VA = shiftdim(sum(py_t(1,1:nRBgoods,:)+py_t(2,1:nRBgoods,:),2))'; %%% total value of RB production
    imports_foreign_RB = vtau.*shiftdim(sum(py_t(3,1:nRBgoods,:),2))'; %% total value of foreign RB sales at home
    RB_import_share = imports_foreign_RB./RB_VA;
    Agg_VA = shiftdim(sum(py_t(1,:,:)+py_t(2,:,:),2))'; %%% total value added
    GDP_2 = shiftdim(sum(py_t(1,:,:)+py_t(2,:,:)));
    GDP_1 = shiftdim(Y_H_t).*shiftdim(P_H_t);
    imports_agg = vtau.*shiftdim(sum(py_t(3,:,:),2))'; %% total value of foreign RB sales at home
    %%Agg_import_share = imports_agg./Agg_VA;
    Agg_import_share = imports_agg./GDP_1';
    %%%plot(year',RB_import_share,'r-',year',0.01*import_shares_RB_wtd_avg_FINAL_data,'k:',year',Agg_import_share,'g-',year',0.01*aggregate_import_share_data ,'k:','LineWidth',2)
    plot(year',RB_import_share,'r-',year',0.01*import_shares_RB_wtd_avg_data,'k:',year',Agg_import_share,'g-',year',0.01*import_shares_US_MFG_data,'k:','LineWidth',2)
    %%axis([1950 2000 0.00 1.00])
    title('Import Shares')
    legend('Rust Belt, Model','Rust Belt, Data','Aggregate, Model','Aggregate, Data','Location','northwest')
    
    %%% Strike rate
    subplot(2,2,4)
    plot(year',shiftdim(100*strike_rate_RB_t)','r-',year',shiftdim(100*strike_rate_SB_t)','b-','LineWidth',2)
    title('Strike Probabilities by Region')
    legend('Rust Belt','Rest of Country')
    ylabel('Percent')
    axis([1950 2000 0 30])
    
    print ('-bestfit','-dpdf','figure_1.pdf')

    
    %%% Figure 2: 3x2 secondary model predictions
    figure;
    hold on;
    subplot(3,2,1)
    plot(year',shiftdim(P_H_t(1,1,:)),'g-',year',shiftdim(P_F_t(1,1,:)),'m--','LineWidth',2)
    title('Price Indices')
    legend('Home','Foreign','Location','northeast')
    
    subplot(3,2,2)
    avg_z_RB = shiftdim(mean(z_t(1,1:nRBgoods,:)))';
    avg_z_ROC = shiftdim(mean(z_t(1,nRBgoods+1:end,:)))';
    plot(year',avg_z_RB,'r-',year',avg_z_ROC,'b-','LineWidth',2)
    title('Average Productivity by Home Region')
    legend('Rust Belt','Rest of Country','Location','northwest')
    
    subplot(3,2,3)
    histogram(log(z_t(1,:,51))-log(z_t(2,:,51)),'Normalization','pdf')
    %%axis([-6 6 0.00 1.00])
    title('Distribution of log zH - log zF in 2000')
    
    subplot(3,2,4)
    avg_z_RB_F = shiftdim(mean(z_t(2,1:nRBgoods,:)))';
    avg_z_ROC_F = shiftdim(mean(z_t(2,nRBgoods+1:end,:)))';
    plot(year',avg_z_RB_F,'r--',year',avg_z_ROC_F,'b--','LineWidth',2)
    title('Average Productivity by Foreign Region')
    legend('Foreign Rust Belt','Foreign Rest of Country','Location','northwest')
    
    subplot(3,2,5)
    prod_growth_RB = z_t(1,1:nRBgoods,2:nPeriods)./z_t(1,1:nRBgoods,1:nPeriods-1);
    prod_growth_ROC = z_t(1,nRBgoods+1:end,2:nPeriods)./z_t(1,nRBgoods+1:end,1:nPeriods-1);
    avg_prod_growth_RB = shiftdim(mean(prod_growth_RB))';
    avg_prod_growth_ROC = shiftdim(mean(prod_growth_ROC))';
    plot(1951:2000,100*(avg_prod_growth_RB-1),'r-',1951:2000,100*(avg_prod_growth_ROC-1),'b-','LineWidth',2)
    %%axis([1950 2000 0.00 4.00])
    title('Average Productivity Growth by Home Region')
    legend('Rust Belt','Rest of Country','Location','southwest')
    
    subplot(3,2,6)
    prod_growth_F_RB = z_t(2,1:nRBgoods,2:nPeriods)./z_t(2,1:nRBgoods,1:nPeriods-1);
    prod_growth_F_ROC = z_t(2,nRBgoods+1:end,2:nPeriods)./z_t(2,nRBgoods+1:end,1:nPeriods-1);
    avg_prod_growth_F_RB = shiftdim(mean(prod_growth_F_RB))';
    avg_prod_growth_F_ROC = shiftdim(mean(prod_growth_F_ROC))';
    plot(1951:2000,100*(avg_prod_growth_F_RB-1),'r-',1951:2000,100*(avg_prod_growth_F_ROC-1),'b-','LineWidth',2)
    title('Average Productivity Growth by Foreign Region')
    legend('Foreign Rust Belt','Foreign Rest of Country','Location','northwest')
    
    print ('-bestfit','-dpdf','figure_2.pdf')
    
else
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% Plot only model
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    figure;
    hold on;
    subplot(3,2,1)
    plot(1:nPeriods,shiftdim(strike_prob_RB_t)','m-',1:nPeriods,shiftdim(strike_prob_SB_t)','b-',1:nPeriods,shiftdim(strike_rate_RB_t)','m--',1:nPeriods,shiftdim(strike_rate_SB_t)','b-')
    title('Strike Probability and Rate')
    legend('Prob, Rust Belt','Prob, Rest of Country','Rate, Rust Belt','Rate, Rest of Country')
    ylabel('Percent')
    axis([1 nPeriods 0 1])
    subplot(3,2,2)
    plot(1:nPeriods,shiftdim(P_H_t(1,1,1:nPeriods)),'m-',1:nPeriods,shiftdim(P_F_t(1,1,1:nPeriods)),'b--')
    title('Price Indices')
    legend('Home','Foreign','Location','southeast')
    subplot(3,2,3)
    plot(1:nPeriods,shiftdim(Y_H_t(1,1,1:nPeriods)),'m-',1:nPeriods,shiftdim(Y_F_t(1,1,1:nPeriods)),'b--')
    title('Output')
    legend('Home','Foreign','Location','southeast')
    subplot(3,2,4)
    plot(1:nPeriods,shiftdim(wF_t(1,1,1:nPeriods)),'m-')
    title('Foreign Wage')
    mean_z_terms = mean(z_t,2);
    mean_z_h_t = shiftdim(mean_z_terms(1,1,:));
    mean_z_f_t = shiftdim(mean_z_terms(2,1,:));
    subplot(3,2,5)
    plot(1:nPeriods,mean_z_h_t,'m-',1:nPeriods,mean_z_f_t,'b--')
    title('Average Productivity')
    legend('Home','Foreign','Location','northwest')
    sum_l_terms = sum(l_t,2);
    labor_h_for_h_t = shiftdim(sum_l_terms(1,1,:))/LS;
    labor_f_for_h_t = shiftdim(sum_l_terms(3,1,:))/LS;
    subplot(3,2,6)
    plot(1:nPeriods,labor_h_for_h_t,'m-',1:nPeriods,labor_f_for_h_t,'b--')
    axis([1 nPeriods 0 1])
    title('Employment Producing for Home Country')
    legend('Home for Home','Foreign for Home','Location','northwest')
    print ('-bestfit','-dpdf','fig_aggregate.pdf')
    
    %% Disaggregate statistics
    
    figure;
    hold on;
    subplot(3,2,1)
    RB_labor_share = shiftdim(sum(l_t(1,1:nRBgoods,:)+l_t(2,1:nRBgoods,:)))';
    plot(1:nPeriods,RB_labor_share,'m-')
    axis([1 nPeriods 0 1])
    title('Rust Belt Employment Share')
    
    subplot(3,2,2)
    avg_z_RB = shiftdim(mean(z_t(1,1:nRBgoods,:)))';
    avg_z_ROC = shiftdim(mean(z_t(1,nRBgoods+1:end,:)))';
    plot(1:nPeriods,avg_z_RB,'m-',1:nPeriods,avg_z_ROC,'b--')
    title('Average Productivity by Home Region')
    legend('Rust Belt','Rest of Country','Location','northwest')
    
    subplot(3,2,3)
    py_t = p_t.*y_t;
    RB_VA = shiftdim(sum(py_t(1,1:nRBgoods,:)+py_t(2,1:nRBgoods,:)))'; %%% total value of RB production
    imports_foreign_RB = shiftdim(sum(py_t(3,1:nRBgoods,:)))'; %% total value of foreign RB sales at homoe
    RB_import_share = imports_foreign_RB./RB_VA;
    plot(1:nPeriods,RB_import_share,'m-')
    title('Rust Belt Import Share')
    
    subplot(3,2,4)
    avg_z_RB_F = shiftdim(mean(z_t(2,1:nRBgoods,:)))';
    avg_z_ROC_F = shiftdim(mean(z_t(2,nRBgoods+1:end,:)))';
    plot(1:nPeriods,avg_z_RB_F,'m-',1:nPeriods,avg_z_ROC_F,'b--')
    title('Average Productivity by Foreign Region')
    legend('Foreign Rust Belt','Foreign Rest of Country','Location','northwest')
    
    subplot(3,2,5)
    avg_i_RB = shiftdim(mean(i_t(1,1:nRBgoods,:)))';
    avg_i_ROC = shiftdim(mean(i_t(1,nRBgoods+1:end,:)))';
    plot(1:nPeriods,avg_i_RB,'m-',1:nPeriods,avg_i_ROC,'b--')
    title('Average Investment by Home Region')
    legend('Rust Belt','Rest of Country','Location','northwest')
    
    subplot(3,2,6)
    avg_i_F_RB = shiftdim(mean(i_t(2,1:nRBgoods,:)))';
    avg_i_F_ROC = shiftdim(mean(i_t(2,nRBgoods+1:end,:)))';
    plot(1:nPeriods,avg_i_F_RB,'m-',1:nPeriods,avg_i_F_ROC,'b--')
    title('Average Investment by Foreign Region')
    legend('Foreign Rust Belt','Foreign Rest of Country','Location','northwest')
    
    print ('-bestfit','-dpdf','fig2_main.pdf')
    
end


